Astronomische Nachrichten, 1 February 2008 
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A simple way to couple an interface dynamo model to a fast tachocline model is presented, under the assumption that 
the dynamo saturation is due to a quadratic process and that the effect of finite shear layer thickness on the dynamo 
wave frequency is analoguous to the effect of finite water depth on surface gravity waves. The model contains one free 
parameter which is fixed by the requirement that a solution should reproduce the helioseismically determined thickness 
of the tachocline. In this case it is found that, in addition to this solution, another steady solution exists, characterized by 
a four times thicker tachocline and 4—5 times weaker magnetic fields. It is tempting to relate the existence of this second 
solution to the occurrence of grand minima in solar activity. 
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1 Introduction 

The cyclically varying magnetic field that gives rise to so- 
lar activity is generally thought to arise as a consequence of 
an a£l dynamo mechanism. The strong toroidal field near 
solar maximum is a product of the solar differential rotation 
that winds up the weak general poloidal field. Poloidal fields 
are then restored by an unspecified "a-effect", the physical 
meaning and spatial location of which has yet to be eluci- 
dated. 

The shear associated with differential rotation is by far 
the strongest in a thin layer below the convective zone, 
known as the tachocline. Current dynamo models therefore 
invariably locate the fi-effect in the tachocline. The site 
and physical nature of the a-effect is much less clear. As 
discussed in recent reviews (Petrovay 2000, Charbonneau 
2005, Solanki et al. 2006), the two most widely duscussed 
alternatives are interface dynamos, where a is concentrated 
near the bottom of the convective zone, separated from the 
tachocline shear region by a thin interface only, and flux 
transport dynamos, where a is concentrated near the sur- 
face, and magnetic flux transport by meridional circulation 
forges a link between the surface and the tachocline. 

Solar dynamo models with a jump in the diffusivity a- 
cross an interface (the bottom of the convective zone) were 
first constructed by Ivanova and Ruzmaikin (1976, 1977). 
Parker (1993) developed an analytical model thas has been 
considered the prototype of interface dynamos ever since. 
Further studies of interface dynamos include Tobias (1996), 
Charbonneau & MacGregor (1997), Markiel & Thomas 
(1999), Petrovay & Kerekes (2005). Flux transport dynamos 
were constructed, among others, by Choudhuri, Schiissler 
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& Dikpati (1995), Dikpati & Charbonneau (1999), Dikpati 
et al. (2004), Chatterjee, Nandy & Choudhuri (2004). 

Like the dynamo models, tachocline models also come 
in two main varieties. Slow tachocline models assume that 
the tachocline lies in the radiative interior, where turbulence 
is insignificant, and it is pervaded by the steady remnant 
magnetic field of the solar interior. If this field is limited to 
the radiative interior, it may be able to confine the tachocline 
to its observed low thickness. Whether or not such a limi- 
tation of the internal magnetic field to the radiative interior 
is feasible, is currently a subject of debate (Sule, Riidiger & 
Arlt 2005, Bran & Zahn 2006). 

Fast tachocline models, in contrast, assume that local in- 
stabilities and/or convective overshoot maintain a moderate 
level of turbulence in the tachocline. In this case the oscil- 
latory magnetic field generated by the dynamo will pene- 
trate the tachocline and determine its dynamics. In fact, as 
current dynamo models place the site of the Sl-effect in the 
tachocline, penetration of the dynamo field in the tachocline 
is a necessary condition for those dynamo models to work. 
Following a suggestion by Gilman (2000), the first detailed 
models of the fast tachocline were constructed by Forgacs- 
Dajka & Petrovay (2001, 2002) and further developed by 
Forgacs-Dajka (2003). 

In fast tachocline models the thickness of the tachocline 
is determined by the intensity of the magnetic field gener- 
ated by the dynamo. The thickness of the tachocline, in turn, 
controls the dynamo via the shear. In this way a nonlin- 
ear coupling is expected to arise between the dynamo and 
a fast tachocline. In the following I will try to explore the 
nature of this nonlinear coupling, for concreteness focusing 
on the case of an interface dynamo. Section 2 describes the 
main features of the simple model used here, while Section 
3 presents the results. Section 4 concludes the paper. 
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2 Elements of the model 
2.1 Interface dynamo 

Parker's analytical kinematic interface dynamo (Parker 
1993) consists of two semiinfinite domains (say, the half- 
spaces z < and z > in an xyz Cartesian frame) of 
diffusivity 77 _ an rj + , respectively. The z > domain is 
characterized by a constant shear (fl = dv y /dz) and no 
a-effect, while the other is characterized by constant a and 
no shear. If the 

~ = ^ »> 

77_?7 + fc J 

dynamo number is sufficiently high, a dynamo wave of form 

(2) 



B oc exp[ikx + (cr + iu})rj + k 2 i\ 



is spontaneously excited at the interface (i.e. a > 0). In the 
limit — > 0, relevant to the Sun, the dimensionless 

growth rate a and the dimensionless frequency w are given 
by the relations 

N = ±8(a+±)(a+l) 1 / 2 t7 1 / 2 (3) 

4a- 1 / 2 if a <C 1 (slightly supercritical case) 
8<7 2 if a 3> 1 (strongly supercritical case) 

and 

a{a + 1) (4) 

a if a <C 1 (slightly supercritical case) 
<t 2 if a ^> 1 (strongly supercritical case) 

In what follows we will focus on the slighly supercrit- 
ical limit. Indeed: interpreting the butterfly diagram as the 
surface manifestation of an interface dynamo wave we find 
for the frequency 



uj = Cb-q+k 2 = 2?r/22years ~ 9 • 1CT 9 s~ 



(5) 



Using A = Rir/2 (R is the radius of the bottom of the con- 
vective zone) we have 

fc = 4/i?~8-10" 11 cm" 1 . (6) 

With a convective zone diffusivity of i] + = 6 • 10 12 cm 2 /s 
(Wang, Nash & Sheeley 1989, Petrovay and van Driel- 
Gesztelyi 1997) this yields Cj ~ 0.23, so a ~ 0.05 < 1. 
Our assumption that the slightly supercritical limit holds in 
the solar case is therefore justified. 

2.2 Finite shear layer 

In Parker's original model, the shear layer representing the 
tachocline is semi-infinite. In reality, the tachocline has 
some finite thickness w. As the difference between the rota- 
tion rate of the radiative interior and the equatorial rotation 
rate of the convective zone is fixed, the mean shear f2 in 
the tachocline, and thereby also N will scale inversely with 
tachocline thickness: N oc 1/w. It is not known how the 
solutions to Parker's problem are influenced if the thickness 
of the shear layer is finite. In analogy to the case of surface 
gravity waves over shallow water here we will assume that 
the finite thickness of the shear layer reduces the squared 



mode frequency uj 2 by a factor tanh(fcw). (At any rate, 
such a factor does have the right asymptotic behavior for 
w — > and w — > oo, expected on physical grounds.) In the 
slightly subcritical case then a is also reduced by the same 
factor, so 



a oc A 2 tanh(fcw) oc tanh(kw) / w 2 



(7) 



In the solar tachocline the angular velocity is expected to 
relax gradually to the rigid rotation rate of the radiative in- 
terior, with some scale height H. Defining w arbitrarily as 
the depth where the residual rotation has been reduced by 
two orders of magnitude, we have w ~ H In 100. Plugging 
this and equation (6) into (7) we have 



a oc tanh(4 In W0H/R) / H 2 



(8) 



Helioseismic studies (Kosovichev 1996, Basu & Antia 
2001) indicate that at low latitude the tachocline is located 
immediately below the adiabatically stratified convective 
zone. The equatorial rotation rate relaxes to the rigid ro- 
tation rate of the radiative interior with a scale height of 
H ~ 5-10 Mm (cf. Forgacs-Dajka & Petrovay 2002); for 
concreteness here we will take H = 7 Mm. This shows that 
the tachocline is very thin indeed, its full thickness w not 
exceeding a few percents of the solar radius. Note that at 
higher latitudes the tachocline may be marginally thicker, 
and it partly overlaps the adiabatic convective zone. 

2.3 Saturation 

Saturation of a kinematic dynamo may be brought about by 
a variety of effects, such as the quenching of the a-effect for 
strong magnetic fields, or magnetic flux loss due to buoy- 
ancy. Without specifying the physical effect responsible for 
it, here we will simply describe the nonlinear saturation in 
a parametric form: 

— =oB- aB 1+h (9) 
at 

This implies that in the saturated state (dB/dt = 0) the 
field strength is 

Bg = a/a = K tanh(41n 100H/R)/H 2 , (10) 

K being an undetermined amplitude. 

For the nonlinearity parameter the most plausible and 
widely used choice is K = 1. (E.g. buoyant flux loss is ex- 
pected to lead to a nonlinear term of this form.) In what 
follows we will therefore focus on the case k = 1. The in- 
fluence of this choice of k on the solution will be discussed 
at the end of Section 3 below. 

Instead of using the heuristic expressions (7)-(10), a 
more satisfactory approach would clearly involve a fully 
consistent solution of the nonlinear dynamo problem with 
a finite tachocline. For certain cases such solutions were 
presented e.g. in Chapter 10 of Zeldovich, Ruzmaikin and 
Sokoloff (1984). However, neither of the cases considered 
there (a-effect and shear located in two infinitesimally thin 
layers or in two semi-infinite layers, respectively, separated 
by a finite distance) does include the effect of a finite layer 
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thickness. An extension of those models to the case of a fi- 
nite shear layer may be a promising way to make further 
advance in the problem, but it is beyond the scope of the 
present paper. 

2.4 Fast tachocline 

The relation between the dynamo field amplitude and the 
tachocline scale height was derived by Forgacs-Dajka & 
Petrovay (2001). Let us recall the main points here. For the 
derivation we consider a plane parallel layer of incompress- 
ible fluid of density p, where the viscosity v and the mag- 
netic diffusivity rj are taken to be constant. At z = a peri- 
odic horizontal shearing flow is imposed in the y direction: 

VyO = v COs(fcx) (11) 

(so that x will correspond to heliographic latitude, while y to 
the longitude). We assume a two-dimensional flow pattern 
(d/dy = 0) and v x = v z = 0. An oscillatory horizontal 
field is prescribed in the x direction as 

B x = B }} cos(u>t) (12) 

The evolution of the azimuthal components of the velocity 
and the magnetic field is then described by the correspond- 
ing components of the equations of motion and induction, 
respectively. Introducing v = v y and using Alfven speed 
units for the magnetic field 

V p = Bp^p)- 1 / 2 b = Byi^p)- 1 / 2 , (13) 

these can be written as 



9v , , db ~ 

- = V p cos(u,t)- + vV*v 



db 



dv 



dt =V * C0S{ut) dx 



i]\7 2 b 



(14) 



(15) 



Solutions may be sought in the form 

v = v(x, z) + v'(x, z)f(u)t) (16) 

b = b'{x,z)f{ut + 4>) (17) 

where / is a 27r-periodic function of zero mean and of am- 
plitude 0(1). (a denotes time average of a, while a' = 

a — a.) 

The (temporal) average of equation (14) reads 

= V p cos(ujt)f(wt + <j>)d x b' - vV 2 v 



Subtracting this from equation (14) yields 



dv 



,db 



— = V p [cos(ujt)f(Ljt + </>)]' — - vV 2 v 



dt 



dx 



(18) 



(19) 



For an estimate, one can suppose cos(wf) / (tut + (/))£ 
0(1) (i.e. no "conspiracy" between the phases, a rather nat- 
ural assumption). As H <C R we may approximate V 2 ~ 
H~ 2 . Estimating the other derivatives as d t ~ to and d x ~ 
R-\ (18) yields 

V p b'/R - vv/H 2 (20) 

A similar order-of-magnitude estimate of the terms in equa- 
tion (19) yields 

(<j + v/H 2 )v' ~ V p b'/R (21) 
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Fig. 1 Tachocline confining field strength B c (dashed) and dy- 
namo saturation field strength B s (solid) as functions of the 
tachocline scale height H for an interface dynamo combined with 
a fast tachocline. Crossing points of the curves represent steady so- 
lutions of the nonlinear dynamo problem. The undetermined am- 
plitude of B 3 has been set so that one solution satisfies currently 
available seismic constraints on H. 



while from (15) we find in a similar manner 

tub' ~ (7J + v')V p /R + r/b'/H 2 



(22) 



From the last three order-of-magnitude relations one can 
work out with some algebra 



V 



_ uR 2 uj (1 + t]/ujH 2 ){1 + u/ujH 2 ) 



(23) 



H 2 1 + 2is/loH 2 

Returning from Alfvenic units to magnetic field strength we 
finally have 

^ ( 1 + ,/ r Hp (1 + : /^) 
H z 1 + 2v/luH z 

For concreteness, in the numerical evaluation of equa- 
tion (24) we will take the following fiducial values: v = 
i] = 10 10 cm 2 /s, p = 0.08 g/cm 3 . 

3 Results 

Setting B = B c = B s , equations (10) and (24) can now 
be solved for H and B. A graphical solution is presented in 
Figure 1 . The dashed line shows relation (24) while the solid 
line represents the relation (10). As this relation involves a 
free amplitude factor K, on this logarithmic plot we are free 
to shift the solid curve in the vertical direction. In the figure, 
the solid curve was shifted so that one crossing point falls 
to the value H = 7 Mm, the most likely value based on 
helioseismic constraints. In this case the mean toroidal field 
strength is seen to exceed one kilogauss (which does not 
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exclude the existence of 10 5 G flux tubes in an intermittent 
field structure, cf. Ruzmaikin 2001). 

It is apparent from Figure 1 that a second stationary so- 
lution also exists, with a roughly 4 to 5 times weaker mean 
magnetic field and a correspondingly thicker tachocline. It 
is tempting to identify this second solution with a "grand 
minimum" state of solar activity, such as the Maunder min- 
imum of the 17th and 18th centuries. From equations (8) 
and (5) we expect that the dynamo period in the grand min- 
imum mode should be about a factor of 2-2.5 longer than 
in normal solar activity. There is indeed some observational 
evidence indicating that during the Maunder minimum the 
dominant period in solar activity may have been around 22 
years (Usoskin & Mursula 2003), in fairly good agreement 
with our expectation. (Note that the present model is plane 
parallel, allowing a continuous mode spectrum. In a finite 
spherical geometry the mode spectrum will be discrete, pos- 
sibly explaining why the dynamo period is exactly doubled 
during grand minima.) 

The results presented above are valid in the case of a 
quadratic nonlinearity, i.e. k = 1. How does a change in the 
value of k influence the validity of the findings? Analytical 
considerations and numerical experimentation show that a 
bimodal solution persists in the range 0.5 < K < 2. How- 
ever, for k < 0.67 the present mode of operation would 
correspond to the weak-field mode of the dynamo, which 
does not seem to agree with historical evidence. On the 
other hand, increasing n above unity, the weak-field solu- 
tion becomes increasingly weaker and its tachocline depth 
diverges as k — > 2. E.g. for k = 1.5 the field strength in 
the weak mode would be two orders of magnitude lower 
than presently, and the tachocline thickness would increase 
to 300 Mm. These values are clearly unrealistic, so the kind 
of bimodal solution outlined above is restricted to satura- 
tion mechanisms whose behaviour is reasonably close to 
quadratic. 

4 Conclusion 

Our model contains one free parameter which is fixed 
by the requirement that a solution should reproduce he- 
lioseismically derived thickness of the tachocline. In this 
case we have found that, in addition to this solution, an- 
other steady solution exists, characterized by a four times 
thicker tachocline and 4-5 times weaker magnetic fields. It 
is tempting to relate the existence this second solution to 
the occurrence of grand minima in solar activity. 

How can the dynamo flipflop from one of these modes to 
the other? Stationary solutions of equation (10) are clearly 
linearly stable for a fixed value of H. Similarly, consid- 
ering the full time-dependent angular momentum equation 
Forgacs-Dajka & Petrovay (2002) found that the fast tacho- 
cline solution is an attractor, i.e. the solution of equation 
(24) for H is also stable if B c is kept fixed. A more gen- 
eral stability analysis of the problem in two variables is yet 
to be done, but we may plausibly expect both solutions to 



be linearly stable if these stationary states are to be realized 
in the Sun for any significant period of time. Switching be- 
tween the two solutions should then be due to finite ampli- 
tude stochastic disturbances. On the other hand, the possi- 
bility of a more complex nonlinear dynamical behaviour of 
the time-dependent system can not be discarded, with deter- 
ministic evolution between two quasi-stationary states. 

The present model clearly hinges on four main assump- 
tions. The dynamo is supposed to operate in an interface 
wave dominated regime (as opposed to a circulation dom- 
inated regime); the saturation mechanism is assumed to be 
quadratic; the effect of finite shear layer thickness on the 
frequency is assumed to be the same as for surface gravity 
waves over shallow water; and finally, the tachocline is as- 
sumed to be confined by the dynamo field. The validity of 
these assumptions still awaits confirmation. 
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